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The  objective  of  this  program  is  to  develop  the  methodology  for  per¬ 
forming  calculations  of  the  turbulent  near  wake  behind  a  slender  body  in 
supersonic  flow.  The  two  major  technical  problems  are: 

1.  To  develop  turbulence  model  equations  which  adequately  describe 
the  mean  flow  properties  and,  in  addition,  provide  information 
about  the  properties  of  the  fluctuation  field. 

2  To  develop  numerical  procedures  for  the  solution  of  the  model 
equations.  f 

A  turbulence  modeling  based  on  an  eddy  viscosity,  which  is  a  function 
of  turbulence  properties  (the  turbulent  kinetic  energy,  e,  and  an  integral 
scale  length  function,  *),  is  being  developed.  The  model  equations  are 
programmed  for  solution  by  a  steady  state  finite  difference  method  (with 
only  the  boundary  layer  form  of  the  laminar  and  turbulent  transport  terms 
included)  for  comparison  with  available  experimental  data.  After  several 
modifications  of  the  original  equations,  the  turbulence  model  has  now  been 
shown  to  adequately  describe  both  mean  and  turbulence  properties  of  incom¬ 
pressible  and  compressible  adiabatic  flat  plate  boundary  loyers.  Additional 
verification  of  the  model  by  comparison  with  mixing  layer  experiments  would 
be  desirable  before  applying  the  model  to  the  near  wake  calculation. 

The  near  wake  problem,  because  of  the  recirculation  region  immediately 
behind  the  body,  constitutes  a  boundary  value  problem.  The  proposed  comput¬ 
ational  approach  is  to  obtain  the  steady  flowfield  as  the  limit  of  the  time 
dependent  solution  of  the  Navier-Stokes  equations  with  turbulence  modeling. 
The  numerical  method  proposed  uses  alternating-direction  implicit  differ¬ 
encing  of  the  conservation  form  of  the  equations.  The  method  has  been 
programmed  for  the  laminar  Navier-Stokes  equations.  Truncation  error  has 
been  shown  to  be  second  order  in  temporal  and  spatial  meshsize.  Stability 
has  been  demonstrated  for  a  temporal  step  size  up  to  ten  times  that  beyond 
which  explicit  difference  methods  become  unstable.  Calculations  of  shock 
structure  have  demonstrated  excellent  agreement  with  exact  solutions  even 
with  relatively  poor  spatial  resolution.  Shock  jump  conditions  are  accur¬ 
ately  calculated  even  with  spatial  resolution  as  poor  as  half  the  shock 
thickness. 
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Aii  initial  appl  ication  or'  the  numerical  methods  to  the  near  wake  pro¬ 
blem,  using  laminar  conservation  equations,  is  in  progress. 
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1 .  INTRODUCTION 


The  principal  objective  of  this  program  is  to  develop  the  capability 
of  performing  calculations  of  turbulent  near  wakes  of  slender  bodies.  Two 
major  subtasks  are  involved: 

1.  Develop  turbulence  model  equations  applicable  to  compressible 
as  well  as  incompressible  flows. 

2.  Develop  finite  difference  methods  for  the  solution  of  the  model 
equations. 

An  integral  part  of  the  development  of  the  turbulence  model  equations 
is  comparison  of  calculations  of  simple  flows  with  experimental  data.  Such 
comparisons  are  necessary  to  test  the  validity  of  the  proposed  modeling. 

These  simple  flows  can  be  modeled  quite  adequately  using  subsets  of  the 
complete  conservation  and  turbulence  model  equations  which  are  being  pro¬ 
posed  for  use  in  the  near  wake  problem.  The  particular  subset  being  used 
for  this  purpose  is  one  which  has  been  applied  to  the  laminar  steady  cal¬ 
culation  of  interacting  flows^,  and  contains  both  the  complete  inviscid 
terms  (in  the  supersonic  flow)  and  the  boundary-layer-1  ike  viscous  terms. 

The  addition  of  the  turbulence  modeling  to  this  existing  code  and  the  develop 
ment  and  testing  of  the  model  equations  using  this  code  has  been  jointly 
supported  by  the  ROPE  Project  funded  by  ABMDA  and  ARPA  under  Contract 
DAHC-60-71-C-0049. 

The  near  wake  problem  has  the  intrinsic  character  of  a  boundary  value 
problem  because  of  the  presence  of  the  subsonic  region  immediately  behind 
the  body.  It  is  possible  to  make  simplifying  assumptions  in  the  recircul¬ 
ation  region  and  recast  the  problem  as  a  pseudo-initial-value-problem.  The 
upstream  influence  characteristic  of  the  boundary-value-problem  is  then 
retained  through  the  presence  of  an  eigenvalue  required  to  pass  through  a 
downstream  saddle-point  singularity  (see  Reference  1).  However,  a  more 
rigorous  approach  (and  the  one  proposed  here)  is  to  describe  the  near  wake 
using  the  complete  steady  state  conservation  equations.  The  boundary  value 
problem  resulting  from  this  formulation  can  be  converted  to  a  more  easily 
solved  Initial  value  problem  by  considering  the  steady  state  solution  to  be 
the  asymptotic  limit  of  a  time-dependent  calculation.  The  numerical  methods 
for  solution  of  these  equations  can  be  developed  relatively  independent  of 


the  turbulence  modeling,  since  the  laminar  (Navier-Stokes)  equations  embody 
most  of  the  computational  difficulties  expected  from  the  complete  set  of 
model  equations. 

This  report  describes  the  initial  (and  relatively  independent)  work 
on  the  modeling  of  turbulent  flows  and  on  the  computational  techniques 
required  to  solve  problems  described  by  Navier-Stokes-1 ike  equations. 

The  work  described  in  this  report  was  performed  in  the  period  1  April 
through  1  October  1972. 
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2.  TURBULENCE  MODELING 


The  calculation  of  high  speed  turbulent  shear  flows  requires  a  descrip¬ 
tion  of  the  Reynolds  stresses  and  the  turbulent  heat  flux  terms  appearing 
in  the  time  averaged  and  hence,  mean  flow  momentum  and  energy  equations. 

Using  arguments  similar  to  those  for  laminar  flows,  the  transport  of  momentum 
and  energy  by  turbulence  is  found  to  be  proportional  to  the  time  averaged 
rate-of-strain  or  time  averaged  temperature  gradient,  and  to  the  so-called 
eddy  diffusivity  function  which  depends  upon  the  structure  of  the  turbulence, 

_ _  /3u.  3U  -  v 

i.e.,  u*|Uj  =  e  ^^7  +  3x7 Classically,  the  eddy  viscosity  has  been 

written  in  terms  of  mean  flow  variables  which  were  found  by  empirical  means 
to  represent  the  turbulence.  For  example,  in  simple  jets,  wakes,  and 
mixing  layers,  e  -  K  6  um=v  where  K  is  an  empirical  constant  and  6  is  the 
lateral  extent  of  the  shearing  region.  The  strong  dependence  of  the  class¬ 
ical  approach  on  the  geometry  of  the  flow  makes  it  unacceptable  for  des¬ 
cribing  a  complex  turbulent  flow  configuration  such  as  encountered  in  the 
near  wake.  Instead,  a  more  sophisticated  representation  of  the  eddy  viscosity 
based  upon  a  local  modeling  of  the  turbulence  field  is  required. 

The  turbulent  field  can  be  represented  by  as  few  as  two  quantities  which 
are  generally  considered  to  be  the  turbulent  kinetic  energy,  e  =  i  u^uT/2 

and  an  integral  scale  length  of  the  turbulence,  £,  which  is  related  to  the 

width  of  the  spectrum  of  the  autocorrelation  function  of  the  turbulent 

kinetic  energy.  The  eddy  diffusivity,  written  in  terms  of  these  variables 
•  1  /  2 

is  e  =  y*  e  £  where  y*  is  a  scaling  constant.  The  deterr'i nation  of  e 
and  «,  requires  the  development  of  two  model  equations  which  describe  the 
convection,  production,  dissipation,  and  diffusion  of  turbulence  in  the 
flow.  An  equation  directly  describing  the  turbulent  kinetic  energy  has 
been  used  in  numerous  modern  investigations,  beginning  with  Bradshaw^, 
is  now  reasonably  well  understood  and  accepted.  However,  less  confidence 
exists  with  the  second  turbulence  model  equation,  primarily  because  this 
equation  represents  a  higher  order  balance  between  production,  dissipation, 
and  diffusion  than  does  the  kinetic  energy  equation.  The  uncertainties, 
which  exist  in  the  modeling  of  individual  terms,  are  therefore  more 


important  in  the  second  equation.  The  determination  of  the  most  appropriate 
second  equation  is  part  of  the  present  effort. 

Two  promising  modeling  theories  for  the  second  turbulence  equation  have 
been  considered  in  the  present  study.  One  was  developed  by  Saffman^  in 
terms  of  a  transport  equation  for  the  turbulent  vorticity  fluctuations,  and 
the  other  was  developed  by  Jones  and  Launder^  who  wrote  a  transport  equa¬ 
tion  for  the  turbulent  energy  dissipation  rate  ed  =  Cd  e3^/£,  where  Cd  is 
a  universal  dissipation  constant  to  be  determined  from  experiment.  To 
determine  which  theory  produces  better  agreement  with  experiment,  both 
equations  were  solved,  each  with  the  turbulent  kinetic  energy  equation,  for 
a  flat  plate  boundary  layer  in  incompressible  flow.  These  equations  are 
mutually  coupled  to  the  mean  flow  momentum  and  continuity  equations  and  the 
resulting  solutions,  therefore,  represent  simultaneous  solutions  of  all 
four  equations. 

The  calculations  were  initiated  by  estimating  on  the  basis  of  experi¬ 
mental  data  the  mean  flow  velocity  profile  and  the  profiles  of  e  and  i 
across  the  boundary  layer.  The  turbulent  boundary  layer  was  then  calculated 
by  streamwise  marching  over  a  distance  of  50  initial  boundary  layer  thick¬ 
nesses.  Experience  of  other  investigators  showed  initial  profile  effects 
to  be  negligible  for  an  integration  of  this  distance.  Referring  to  the 
vorticity  model  of  Saffman,  by  adjusting  the  constants  appearing  in  the 
turbulence  equations  over  their  acceptable  range,  it  was  possible  to  produce 
solutions  which  were  in  reasonable  agreement  with  experimental  data  in  the 
law-of-the-wal 1  region.  Two  such  cases  are  shown  in  Figure  1  where  the 
velocity  ratio  U/U  (U  =  (t  /p  )ly^2)  is  plotted  against  the  reduced  coordin- 
ate,  y*  =  p  u  y/u  .  A  wealth  of  experimental  data  has  been  correlated  in 

W  T  W 

the  law-of-the-wal 1  region  with  the  expression 

U_  =  C  +  1  ^n  y* 

T 

and  the  vorticity  model  results  are  in  reasonable  agreement  with  this  cor¬ 
relation  with  C  =  5.24  and  <  =  0.40  in  the  range  100  <  y*  <  2000.  However, 

(5) 

experimental  data  by  IJeighardt'  is  also  presented  and  shows  a  breakaway 
from  the  law-of-the-wall  region  at  about  y*  =  800,  indicative  of  the  outer 
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wall-of-the-wake  region.  The  vorticity  model  thus  fails  to  produce  the 
wal 1-of-the-wake  region  of  the  boundary  layer  consistent  with  experiment. 

This  is  more  clearly  seen  in  Figure  2  where  a  comparison  of  one  vorticity 
model  solution  in  velocity  defect  coordinates  can  be  made  with  the  experi¬ 
mentally  correlated  velocity  defect  law.  The  vorticity  model  solution  does 
not  display  the  curvature  required  to  agree  with  the  velocity  defect  law  in 
the  wake  region  (y/6  >  0.2). 

Further  examination  of  the  solution  shows  in  Figures  3  and  4  that  the 
integral  scale  length,  a/6,  and  dimensionless  kinematic  eddy  viscosity 
e/ue6*,  where  6*  is  the  displacement  thickness,  both  reach  values  in  the 
central  portion  of  the  boundary  layer  which  far  exceed  the  measurements  of 
Klebanoff^  and  the  values  deduced  by  Maise  and  McDonald^  from  experi¬ 
mentally  measured  and  correlated  velocity  profiles.  A  comparison  of  the 
calculated  turbulent  kinetic  energy  is  made  with  Klebanoff's  data  in  Figure 
5,  and  the  agreement  is  clearly  much  better. 

Since  £  ^  e^A,  the  discrepancies  noted  in  Figure  4  between  the 
theoretical  (vorticity  model)  and  measured  e  are  mostly  attributable  to 
the  length  scale,  a.  This  suggests  that  the  problem  lies  with  the  vorticity 
equation.  Physically,  the  difficulty  seems  to  be  attributable  to  the  inabil¬ 
ity  of  the  turbulent  vorticity  fluctuation  function  to  diffuse  from  the  outer 
wake  region  towards  the  wall;  consequently,  the  wall  has  much  too  large  an 
influence  in  the  central  portions  of  the  boundary  layer. 

This  calculation  was  repeated  with  the  energy  dissipation  rate  model 
equation  of  Jones  and  Launder  and  the  results  were  founr^  to  be  in  accept¬ 
able  agreement  with  experimental  data.  The  velocity  distribution  U/U^ 
versus  y*  is  presented  in  Figure  1  and  clearly  displays  law-of-the-wall 
and  law-of-the-wake  regions  which  are  in  good  agreement  with  Weighardt's 
data,  and  the  correlation  curve.  Some  indication  of  the  sensitivity  of 
the  solution  to  the  constants  is  also  given.  Curves  C  and  D  differ  oecause 

of  a  small  change  in  the  production  constant  C_  and  sublayer  constant  C  . 

9  b 

In  Figure  2,  the  agreement  of  the  velocity  profile  from  the  dissipation 
model  with  the  velocity  defect  law  is  excellent,  and  demonstrates  the  ability 
of  the  theory  to  predict  the  law-of-the-wake  region.  The  integral  scale 
length,  A/6,  the  kinematic  eddy  viscosity,  e/Ue5*,  and  the  turbulent  kinetic 


Figure  3.  Comparisons  with  Experiment  of  Calculated  Scale  Length 
Profiles  for  Mach  Zero  Adiabatic  Flat  Plate  Flow 


Figure  4.  Comparisons  with  Experiment  of  Calculated  Kinematic  Eddy 
Viscosity  Profiles  for  Mach  Zero  Adiabatic  Flat  Plate  Flow 


energy  function,  /2e/UT ,  In  Figures  3,  4,  and  5,  respectively,  are  also  In 
much  better  agreement  with  the  data,  and  the  differences  which  are  present 
are  not  considered  serious. 

Jones  and  Launder  applied  thle  dissipation  model  to  accelerating  tur¬ 
bulent  boundary  layers  undergoing  relamlnarizatlon  and  obtained  results 
which  were  in  remarkably  good  agreement  with  existing  data.  Since  an 
Important  part  of  the  present  task  deals  with  an  accelerating  boundary  layer 
about  the  rounded  aft  shoulder  of  a  reentry  vehicle  prior  to  separation, 
their  results  are  very  encouraging.  However,  an  extremely  Important  facet 
of  the  modeling  Is  the  Inclusion  of  the  p>*oper  compressibility  effects,  and 
hence  the  next  step  In  the  modeling  was  In  this  direction. 

Two  experimentally  observed  features  of  compressible  turbulent  flat 
plate  flow,  with  which  the  theoretical  model  must  agree,  are  that  the  velo- 

r  /0  \]/2 

city  profile  under  the  transformation  u*  =  J  I  du  and  the  integral 

o  e 

length,  i,  are  essentially  Invariant  with  Mach  number.  An  analysis  of  the 

law-of-the-wall  region  for  compressible  flow  demonstrated  that  the  two 

3/2 

turbulence  model  equations  written  for  the  variables  pe  and  p  ed  should 
yield  the  desired  Invarlence  of  i  with  Mach  number.  This  was  confirmed 
in  a  calculation  of  a  Mach  5  adiabatic  boundary  layer,  and  the  results  are 
shown  <n  Figure  6  where  a  comparison  with  the  Mach  zero  (Incompressible 
case)  Is  made.  The  results  of  two  other  attempts  where  the  variables  were 
e,  p3/2cd,  and  e,  ed/p  are  also  shown,  and  are  In  strong  disagreement  with 
the  f'jch  zero  results.  A  comparison  of  the  turbulent  kinetic  energy  with 
experimental  data  of  Klstler^  Is  presented  In  Figure  7  and  shows  that 

3/5 

acceptable  agreement  with  the  variables  pe  and  p  ed  is  achieved.  Better 
agreement  Is  obtained  with  the  variables  e  and  ed/p,  but  the  large  discrep¬ 
ancy  In  i  prohibits  use  of  this  function.  The  effect  of  Mach  number  on 
the  eddy  viscosity  Is  Illustrated  In  Figure  8,  both  from  the  calajpatlon 
and  from  the  experimentally  deduced  results  of  Malse  and  McDonald.  The 
predicted  trend  with  Mach  number  Is  correct  even  though  the  values  are  some¬ 
what  larger  than  the  Malse  and  McDonald  results.  A  comparison  of  the  velo¬ 
city  profiles  at  Mach  5  and  Mach  zero  Is  shown  In  Figure  9  under  the 
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Figure  5.  Comparisons  with  Experiment  of  Calculated  Turbulent  Kinetic 
Energy  Profiles  for  flach  Zero  Adlalwtic  Flat  Plate  Plow 
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Figure  6.  Turbulent  Scale  Lengths  for  Different  Compressibility  Transformations 
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Figure  8.  Effect  of  Mach  Number  on  Eddy  Viscosity  Dissipation  Model  <?e 
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Figure  9.  Turbulent  Modeling  Comparison  of  Mach  5  and  Mach  Zero  Velocity  Profiles  Adiabatic  Flat  Plate 


aforementioned  transformation.  The  two  curves  should  be  In  better  agreement; 
however,  the  early  breakaway  of  the  Mach  5  profile  from  the  laminar  sublayer 
(y*«  20)  Is  due  to  a  compressibility  effect  not  accounted  for  in  the  sub¬ 
layer  formulation.  The  turbulent  viscosity  used  In  the  results  shown  is 
modified  by  the  local  turbulent  Reynolds  number  Re  =  pe^t/u  according  to 

uT  =  pe1/2*/(l  +  .224  Cs/Re£) 

Subsequent  analysis  has  shown  that  the  factor  (p/pw)2  should  multiply  the 
constant  Cs,  and  that  the  effect  should  be  to  bring  the  velocity  profile 
Into  good  agreement  with  the  Mach  zero  profile.  A  recalculation  of  the 
boundary  layer  with  this  minor  change  Is  anticipated  in  the  immediate 
future. 

The  development  of  the  turbulence  modeling  has  taken  several  twists 
and  turns,  and  further  comparisons  with  experiment  are  yet  required.  For 
example,  how  the  modeling  holds  up  under  a  mixing  layer  calculation  is 
Important  to  the  near  wake  problem.  The  decay  of  the  turbulent  remnant  of 
the  boundary  layer  Is  also  Important  during  the  expansion  into  the  near 
wake,  and  It  seems  likely  that  additional  modifications  will  be  made  to 
Improve  the  model  for  pressure  gradient  effects.  The  turbulence  model 
equations  as  thpy  stand  today  are  presented  below  (minus  pressure  gradient 
effects).  The  dissipation  rate  equation,  as  given,  is  the  one  used  to 
obtain  the  results  presented  here. 


Turbulent  Kinetic  Energy  Equation 


(conv.)  (prod.) 

(  3u\2 


1  A. 


(diss. ) 
£2  1 


\ 


IY 


(diff.) 

\  . 


1 


(conv. ) 


(prod. ) 


(diss. ) 


where  n:  distance  normal  to  wall 

yT  =  pe1/2a/(l  +  .224  Cs/R£) 

ed  =  Cd  e3/2/£ 

R£  =  pe1/2£/y 

and  the  constants  are: 


Cd  =  0.09 

C  =  120 
s 

Cg  =  1.090 

Q 

II 

<s*4 

Cm  =  0.055 
m 

o*  =  1.0 

Ci  =  2.5  -  Cg 

C2  =  2.5  - 
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3.  NUMERICAL  METHODS 


The  TACIT  (Temporal  Alternating-direction-implicit  code  for  Compressibl 
and  Incompressible  Turbulent  flows)  code  is  being  developed  to  solve  the 
complete  turbulence  model  equations.  As  a  first  step,  the  code  has  been 
written  to  solve  the  compressible  time-dependent  Navier-Stokes  (laminar) 
equations.  This  section  describes  and  examines  the  characteristics  of  the 
finite  difference  methods  used. 

A  number  of  numerical  techniques^  ^have  been  applied  to  related 
problems.  With  the  exception  of  one  method^1  3^,  these  use  explicit  differ¬ 
ence  approximations,  which  are  relatively  easy  to  apply,  but  which  have  in 
common  a  stability  requirement  of  the  form 


where 


max  [aKc,  b«R]  <  1 

(1) 

„  _  (u  +  c) 

Kc  •  —&(—  4t 

(2) 

„  _  yAt 

R 

(3) 

a  and  b  are  constants  of  order  unity 
c  =  local  speed  of  sound. 

This  limitation  on  the  time  step  size  makes  these  methods  impractical  for 
lnw  speed  flows  (u  «  c),  where  the  time  interval  must  be  small  enough  that 
a  sound  wave  travels  less  than  one  spatial  mesh  interval.  For  high  speed 
flows,  this  limitation  is  not  generally  as  serious,  except  that  high  spat¬ 
ial  resolution  (small  AX)  becomes  costly  since  At  must  then  also  be  reduced. 
This  can  become  a  serious  limitation  when  only  limited  regions  of  a  flow- 
field  require  a  small  spatial  meshsize.  This  is  particularly  true  of 
turbulent  boundary  layers,  which  require  very  high  resolution  for  the  lam¬ 
inar  sublayer. 

The  one  partially  implicit  scheme^13)  which  has  been  applied  to  related 
problems  is  better  s ji ted  than  the  explicit  methods  to  low  speed  flows,  but 


has  stability  requirements  which  make  it  unsuitable  for  describing  viscous 
effects  at  high  Reynolds  numbers 

A  method,  which  has  not  been  applied  to  the  Navier-Stokes  equation, 
but  appears  to  be  well  suited  for  this  application  for  both  high  and  low 
speed  flows,  with  no  stability  limitations  and  .asonably  small  truncation 
errors,  is  the  alternating-direction-implicit  (ADI)  scheme.  This  method 
was  first  applied  to  the  solution  of  the  two-dimensional  heat  conduction 
equation 


3U  92u  . 
3t  3? 


S2U 


(4) 


by  Peaceman  and  Rachford^17^  and  to  a  system  of  hyperbolic  equations  of  the 
form 


3u  R 

al+  B 


9U 

sy 


(5) 


(18) 

by  Gourlay  and  Mitchelr  ,  and  was  shown  to  be  unconditionally  stable  in 
both  cases.  The  method  would  be  expected  to  have  the  same  stability  char¬ 
acteristics  when  applied  to  the  Navier-Stokes  equations.  The  following 
sections  briefly  describe  the  ADI  difference  approximations  and  their 
application  to  the  present  problem.  A  more  detailed  description  will  be 
presented  in  a  forthcoming  special  report. 

3.1  Equations 

The  TACIT  code  is  designed  to  solve  the  continuum  conservation  equations 
in  rectilinear  axisymmetric  coordinates.  These  equations  can  be  written 
in  the  form: 

continuity 

H + f!f  (rpv)  =  0  (6) 
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axial  momentum 


i£ii  +  i£H_+  1  i_  (rpUV)  +  If- 
3t  az  r  3r  ^rpuv;  az 


3o 


ii+  li_  (re  )  +  f 
az  r  ar  v  rz  z 


radial  momentum 


M.  +  i£!il  +  1  i_  (rPv2)  ^  = 

at  az  r  ar  upv  ;  ar 


90 rz  +  1_  a_  /  v  fee 
az  r  ar  ^rorr' 


r  +Fr 


total  energy 

M+ ifiHH+li/rpHv)  =  i_[.q  +  Uo  +  vo  ] 
at  az  r  ar  v  p  '  gz  1  Hz  zz  rzJ 

+  —  —  fr(-q  +  uo  +  Vo  )]  +  uF  +  vF„ 
r  sr  L  '  t  rz  rr,J  z  r 

All  quantities  are  non-dimensional ,  using  the  following  reference 

F  =  FL/pruJ 

u2  v2  _  2 

E  =  e  +  ^  +  2  =  F/u^ 

2  2 

H  =  h  +  ^  +  7T  =  fl/u^  where  h  =  e  +  - 
2  2  r  p 


q  =  q/prur 


z  =  z/L 


r  =  r/L 


t  =  t  ur/L 
p  =  p/pr 
u  =  u/ur 
V  =  v/ur 


C7J 


(8) 


(9) 


values : 
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P  =  F/p  u? 

r  r  1 

.  :  I  |  I 

-z  2  i 

a  =  o/Prur  , 

,  i 

The  stresses  a..  act  on  the  normal  p^ne  i  in  the  direction  j.  The  heat1 
i  J 

flux  q.  is  in'  the  direction  i.  The  bodv  force  F.  is  in  the  direction  i. 
The  energy  equation  is  not  completely  general1,  in.  that  radiation  is  not 
included;  mul ti -component  gases  (with  or  without  chemical  reactions)  are 
properly  described  only  if  all  diffusion  poefficients  are  equal  and  the 
Lewis  number  is  unity,  and  if  in  addition  diffusion  is  due  to  gradients  in 
species  concentration  only. 


As  an  initial  test  of  the  computing  technioue,  the  conservation 
equations  are  programmed  for  solution  of  £  laminar  flow,  v/ith  no  body 
forces  present,  so  that: 
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j(2u 

+  \) 

JS\ 

DZ 

X  D_ 

r  3r 
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3r) 

1 

i 

_n  =  _J±_  JiL 
Mz  RePr  ;3z 


_n  =  _J L_  ill 

Hr  RePr  3r 


(10) 

(11) 

(12) 

(13) 

(14) 
(1  5) 


Finally,  it  is  assumed  that  an  equation  of  state  is  available,  of  the 
form: 


(P,h) 


(16) 
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3.2  Difference  Approximations 

Two  step  difference  approximations  are  used,  with 

(g  •  (**'  -  Ak  V 

\3t'1,j  \i,j  i.j/ 


(17) 


(•si, 

t'&L  '  Bi.j  'A,-J-’)/24r 


(18) 


(19) 


J 


j+l  '  2Ai , j  +  ?  A1 


Y 

.j-i/ 


/Ar 


(20) 


(21) 


(»  Cu = Vi,j  (•  + 2a'-.j  •  ?Ai-i.j)k+1/“ 

(B  "L-i  *  Vi  (•?Au*'t2A'.J^A*.J-')k/4r 


(22) 


(23) 


for  the  first  (odd)  time  step.  Here, 


tk+l  .  tk  + 


At 


ziti  !  VAZ 


rjtl  ,ri!lr 
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The  even  time  step  approximations  are 


(24) 


(25) 


(26) 


Here,  for  the  odd  step,  derivatives  in  the  z  direction  are  implicit  (written 
in  terms  of  dependent  variables  at  the  new  time  point)  and  derivatives  in 
the  r  direction  are  explicit.  For  the  even  step,  z  derivatives  are  explicit 
while  r  derivatives  are  implicit.  When  these  difference  approximations  are 
applied  to  the  Navier-Stokes  equations  in  conservation  form,  the  result  is 
a  set  of  algebraic  non-linear  equations  for  values  of  the  dependent  variables 
at  time  k+1.  Coupling  between  values  at  1+1 ,  i,  1-1  occurs  during  the  odd 
step  so  that  equations  for  all  values  of  1  are  coupled  together  for  any 
particular  value  of  j,  and  must,  therefore,  be  solved  simultaneously. 
Conversely,  during  the  even  step,  the  equations  for  all  values  of  j  must 
be  solved  simultaneously  for  any  particular  value  of  i.  The  coupled  non¬ 
linear  algebraic  equations  are  solved  by  a  quasi-linearization  technique. 
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3.3  Truncation  Error 


t;j 


j 


Analysis  shows  that  the  Navier-Stokes  equations,  differenced  using  the 
approximations  of  Section  3.2,  have  a  truncation  error  of  order  At2,  az2, 
or  .  This  has  been  verified  for  the  Invlscld  terms  by  performing  a  calcula¬ 
tion  In  which  a  spherically  synmetrlc,  Gaussian  shaped  energy  pulse  Is 
Instantaneously  Introduced  Into  an  otherwise  uniform  very  large  Reynolds 
number  flow.  An  analytic  solution  can  he  obtained  for  this  problem  and  the 
calculated  error  Is  Indeed  second  order  In  at,  Ar»  Ar. 

While  It  Is  clear  that  a  small  truncation  error  Is  preferable  over  a 
large  one,  It  Is  usually  not  obvious  how  large  a  truncation  error  can  be 
tolerated  before  the  solution  Is  unacceptably  degraded.  In  fact,  very  large 
truncation  errors  may  be  acceptable  for  problems  In  which  error  accumulation 
Is  not  of  overriding  Importance  (such  as  boundary  layer  problems,  where 
boundary  conditions  dominate  the  solution  and  the  history  of  the  flow, 
Including  previously  generated  errors,  becomes  progressively  less  important 
with  distance  downstream).  On  the  other  hand,  wave  Interaction  problems, 

In  which  the  phase  of  a  wave  (which  propogptes  over  many  wavelengths  within 
the  domain)  Is  Important,  require  extremely  small  truncation  errors. 

In  order  to  crudely  simulate  the  requirements  of  the  near  wake  problem, 
the  following  problem  was  numerically  solved: 


t  ■  0:  h  -  1  +  .005  e 


[(z-z  )2  ♦  (r-r  )2]/4 


P 


p  •  1,  u  ■  1,  v  *  0 


where  P  »  ^  ph  (Ideal  gas) 

M  E  4,  Re  ■  2 COO 

For  large  rQ,  this  corresponds  to  an  Instantaneous  cylindrical  Gaussian 
shaped  energy  pulse  at  time  zero  In  a  Mach  4,  otherwise  uniform  flow.  The 
translational  motion  can  be  transformed  out  so  the  resulting  disturbance 
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Is  one-dimensional  (cylindrical),  centered  about  the  point  r  ■  rQ,  2  *  zQ  ♦  t. 
Contours  of  scalar  quantities  (Isotherms,  Isobars,  etc.)  are  therefore  cir¬ 
cular  In  shape,  and  the  degree  of  distortion  of  these  contours  Is  a  quick 
visual  measure  of  the  accumulated  errors  In  the  calculation.  It  was  found 
that  at  ■  0.8,  az  ■  Ar  ■  1  produced  a  solution  which  Is  marginally  accept¬ 
able  after  16  time  steps  (t  -  12.8).  The  degradation  of  this  solution  with 
Increasing  At  or  Increasing  az  and  Ar  Is  Illustrated  In  Figure  10,  which 
shows  contours  of  |  P  -  -J-y  ),  the  pressure  disturbance.  The  most  serious 

\  J 

manifestation  of  accumulating  truncation  error  Is  the  appearance  of  a  wake 
of  error  trailing  behind  the  disturbance.  This  wake  consists  of  cells  of 
alternating  sign,  with  the  largest  error  associated  with  the  cell  closest 
to  the  disturbance.  One  cel i  Is  added  with  each  time  step,  so  the  plot  of 
the  solution  with  the  fewest  time  steps  (At  ■  6.4,  two  time  steps)  has  an 
uncluttered  appearance,  but  with  large  errors. 

The  conditions  of  this  problem  are  In  some  ways  more  severe  than  would 
be  encountered  In  the  near  wake  problem.  The  spatial  resolution  of  the 
Gaussian  disturbance  Is  quite  poor  and  these  conditions  would  be  paralleled 
In  the  near  wake  problem  only  near  the  body  corner,  where  a  strong  expansion 
produces  large  gradients,  and  at  the  separation  and  wake  shocks.  However, 
these  regions  are  stationary  or  at  most  slowly  varying  with  time  and  can, 
therefore,  be  expected  to  have  much  smaller  temporal  contributions  to  the 
error. 

It  appears  that  doubling  At  or  doubling  az  and  Ar  over  the  values  at 
the  nominal  conditions  (At  *  0.8,  ax  *  Ar  ■  1)  have  a  similar  effect  on  the 
magnitude  of  the  errors.  It  Is  reasonable  to  assume,  therefore,  that  the 
spatial  and  temporal  contributions  to  truncation  error  are  of  the  same  order 
of  magnitude  under  these  conditions.  For  the  particular  problem  being 
solved,  one  should  expect  little  further  Improvement  In  the  solution  If  only 
At  Is  decreased,  leaving  ax  and  Ar  fixed.  On  the  other  hand,  for  a  problem 
which  Is  approaching  a  stationary  state,  the  temporal  contribution  to  trunc¬ 
ation  error  will  decrease  considerably,  and  the  time  step  size  can  be  corres¬ 
pondingly  Increased  without  a  significant  degradation  of  the  solution. 

Although  It  Is  unrelated  to  truncation  error,  It  should  be  noted  here 
that  the  nominal  conditions  In  Figure  10  (At  ■  0.8,  ax  -  Ar  »  1)  correspond 
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to  a  Courant  nunfcer  Kc  ■  1.  With  Increasing  Ac,  Kc  ranges  from  1  to  8, 
with  no  sign  of  Instability.  It  therefore  appears  that  In  the  near  wake 
calculation,  especially  In  the  later  stages  as  the  solution  approaches  a 
stationary  state,  a  Courant  number  considerably  greater  than  unity  can  be 
used  without  serious  degradation  of  the  solution  due  to  truncation  error. 

3.4  Calajla tlon  of  Shocks 

The  differencing  of  the  Invlscid  terms  In  their  conservative  form  is 
equivalent  to  an  Integration  of  the  differential  equation.  The  Integrated 
conservation  equations  are,  therefore,  satisfied  even  if  the  detailed  pro¬ 
files  between  two  me'sipolnts  are  unresolved.  This  Is  most  Important  when 
Internal  shocks  appear  In  the  domain  of  the  calculation,  since  the  shock 
jump  conditions  are  then  accurately  calculated.  This  property  of  the  method, 
as  well  as  stability  and  accuracy  charicterlstlcs,  can  be  tested  by  compar¬ 
ing  the  numerical  solution  of  a  norma'  shock  layer  problem  with  the  exact 
solution  obtained  by  G.  1.  Taylor  for  the  limiting  case  of  Infinite  Prandtl 
number. 

In  this  limiting  case,  the  one-dimensional  equations  In  conservation 
form  are  given  vectorlally  by 

$Hf‘°  «27» 

pU 

,  g  -  P  +  PU2  -  ozz 
puH  -  uoz2 

with 

E  -  H  -  P/p 
•zz  ■  fe  <2“  ♦  »  31 

) 


27 


Taylor's  exact  steady  state  solution  of  Equations  27,  for  an  Ideal  gas 
equation  of  state,  can  be  written 

(a  +  ku)A  (b  -  ku)  -  eBz  (28) 


where 


*-(#)[' 

■■(' -*)(>*) 


’1 


2  +  (y-1 )M.‘ 

- - g— 


1  -  M 

I 

(y+Dm/ 

(y+1 

^(M^-l) 


Uj  and  are  respectively  the  velocity  and  Mach  number  ahead  of  the  shock, 
while  m  Is  the  mass  flow  rate  across  the  shock,  pu.  Equation  28  has  been 

U1  +  u2 

normalized  In  such  a  way  that  u  ■  — ^ —  at  z  *  0. 

The  one-dimensional  version  of  the  differencing  presented  In  Section 
3.2  applied  to  Equations  27  gives  difference  equations  of  the  form 


odd  cycle 


1 


«[ 

even  cycle 

1 


k+1 


1  i  r  k+1  k+1  l 

-  f  +  577  9  -  9  =  0 

J  zaz  L  1+1  1-1  J 


At 


r k+i  i 

f  -  f  +  577  9  -  9 

l  J  zaz  L  1+1  1-1  J 


(29) 


(30) 


where  the  Indices 


1,  k  are  Implied  unless  otherwise  specified. 


We  observe  that  the  odd  cycle  Is  implicit  while  the  even  cycle  is 
explicit.  The  solution  of  Equation  29  thus  requires  the  quasi-linearization 
of  the  non-linear  components  of  the  functions  f  and  g;  and  the  solution  by 
Iteration  of  a  matrix  equation  of  the  type 

flw  ■  B" 

where  fl  is  a  band  matrix  with  a  width  of  11  elements. 

The  resolution  and  Courant  number  effects  on  the  stability  and  accuracy 
of  the  numerical  scheme  were  Investigated  by  solving  Equations  29  and  30  to 
steady  state  and  comparing  with  Taylor's  exact  solution. 

Some  results  of  calculations  of  shock  structure  for  (the  Mach  number 
of  the  undisturbed  flow  relative  to  the  shock)  =  2  are  shown  In  Figures  11 
and  12.  In  Figure  11,  the  Courant  number  was  fixed  at  0.75.  The  three 
numerical  solutions  shown  Illustrate  the  effect  of  spatial  resolution.  The 
shock  thickness  for  these  conditions  corresponds  to  d-j 6/vj-j  k10.  The 
solutions  have  been  shifted  In  z  to  give  a  best  fit  with  the  exact  solution 
since  the  final  stationary  position  of  the  shock  within  the  calculation 
domain  Is  not  relevant  In  making  these  comparisons.  The  calculation  with 
a  spatial  meshslze  p^^az/p^  =  5.23  therefore  leaves  the  shock  structure 
entirely  unresolved.  The  jump  conditions  in  velocity  (and  those  In  the 
other  flow  variables  which  are  not  presented  here)  are  nevertheless  accur¬ 
ately  computed.  The  shock  structure  Is  already  represented  reasonably  well 
when  the  mesh  spacing  Is  halved  to  2.62.  The  undershoot  in  velocity  at  the 
downstream  edge  of  the  shock,  which  Is  often  observed  in  other  finite 
difference  methods  and  can  be  quite  large.  Is  not  significant  here. 

In  Figure  12,  the  resolution  has  been  fixed  at  p^u^az/p^  =  1.31.  The 
time  stepslze  has  been  varied  to  examine  the  effect  of  Courant  number  on 
stability.  The  calculation  was  found  to  be  stable  even  for  a  Courant  number 
as  large  as  10.  The  temporal  truncation  error  associated  with  this  large 
time  step  size  apparently  does  not  prevent  the  approach  to  a  steady  solution. 
Since  It  Is  likely  that  an  oscillatory  or  divergent  solution  which  fails  to 
to  approach  a  steady  state  would  result  when  the  time  step  size  becomes 
large  enough  that  the  temporal  truncation  error  dominates  the  solution,  It 
can  be  assumed  that  the  temporal  truncation  error  Is  still  reasonably  small 
even  with  this  large  step  size. 
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*CIT  CODE  FINITE  DIFFERENCE  SOLUTIONS 


Figure  11.  Effect  of  Spatial  Resolution  on  Computed  Shock  Structure 


Figure  12.  Effect  of  Temporal  Step  Size  on  Shock  Structure  Calculation 


Since  explicit  difference  methods  are  limited  by  stability  considera¬ 
tions  to  K  <,  1 »  the  calculation  for  K£  =  10  reaches  a  steady  state  in 
approximately  one  tenth  the  number  of  time  steps  required  by  an  explicit 
method,  with  no  significant  degradation  of  the  steady  solution  due  to  tem¬ 
poral  truncation  error.  This,  again,  underlines  the  fundamental  advantage 
of  the  implicit  method,  especially  when  the  time  dependent  calculation  is 
being  used  to  obtain  the  asymptotic  steady  solution. 


4.  LAMINAR  NEAR  WAKE  CALCULATION 


A  high  Reynolds  number  slender  body  near  wake  calculation  is  under 

progress.  The  body  is  an  8  degree  half-angle  sharp  cone,  with  =  21, 

Re  =  7.8  x  106,  H, ,/H  =0.02  (cold  wall).  Under  these  conditions,  the 
“p  ’  w  <» 

boundary  layer  at  the  rear  of  the  body  is  expected  to  be  fully  turbulent. 

A  steady  calculation,  using  the  methods  of  Reference  1,  has  previously  been 
made  for  these  conditions  under  the  assumption  that  while  the  boundary 
layer  profiles  at  the  body  shoulder  correspond  to  boundary  layer  transition 
near  the  nose,  the  mixing  in  the  near  wake  is  dominated  by  the  laminar 
viscosity.  For  the  purpose  of  further  testing  the  numerical  techniques 
for  the  time- dependent  calculation  under  conditions  very  close  to  those 
correspondin'-  +o  the  final  (turbulent  mixing)  case,  this  calculation  is 
being  used  as  a  test  case  for  the  laminar  version  of  the  code. 
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